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Abstract 

We present the first high order one-step ADER-WENO finite volume scheme with 
Adaptive Mesh Refinement (AMR) in multiple space dimensions. High order spatial 
accuracy is obtained through a WENO reconstruction, while a high order one-step 
time discretization is achieved using a local space-time discontinuous Galerkin pre- 
dictor method. Due to the one-step nature of the underlying scheme, the resulting 
algorithm is particularly well suited for an AMR strategy on space-time adaptive 
meshes, i.e. with time-accurate local time stepping. The AMR property has been 
implemented 'cell- by-cell', with a standard tree-type algorithm, while the scheme 
has been parallelized via the Message Passing Interface (MPI) paradigm. The new 
scheme has been tested over a wide range of examples for nonlinear systems of hy- 
perbolic conservation laws, including the classical Euler equations of compressible 
gas dynamics and the equations of magnetohydrodynamics (MHD). High order in 
space and time have been confirmed via a numerical convergence study and a de- 
tailed analysis of the computational speed-up with respect to highly refined uniform 
meshes is also presented. We also show test problems where the presented high order 
AMR scheme behaves clearly better than traditional second order AMR methods. 
The proposed scheme that combines for the first time high order ADER methods 
with space-time adaptive grids in two and three space dimensions is likely to be- 
come a useful tool in several fields of computational physics, applied mathematics 
and mechanics. 

Key words: Adaptive Mesh Refinement (AMR), time accurate local timestepping, 
space-time adaptive grids, High order WENO reconstruction, ADER approach, 
local space-time DG predictor, hyperbolic conservation laws, Euler equations, 
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1 Introduction 



The idea of using fully discrete one-step time update schemes, rather than 
multi-stage Runge-Kutta time integrators, dates back to van Leer [99j and 
Harten et al. [46J, who applied it to the case of MUSCL and ENO schemes, 
respectively. The key feature of one-step time update schemes is that, after 
a piecewise polynomial approximation to the solution has been obtained at 
time t n through a prescribed reconstruction procedure, an element-local time 
evolution of such reconstructed polynomials is performed, thus allowing for a 
better than first order computation of the numerical fluxes between adjacent 
cells. The one-step time integrators can also be seen as a particular procedure 
to solve (approximately) the generalized Riemann problem at the element in- 
terfaces, where the initial data consists in piecewise polynomials instead of 
piecewise constant data, as it was the case in the original first order Go- 
dunov scheme [Hj. The generalized Riemann problem has been discussed in 
[38TT 6 , 39 . 95 1 191165] and has been used as a building block of numerical methods 
e.g. in [91TL0] as well as in the ADER approach [92,96,90,87,35,32,29,3,5] . In 
the original version of the ADER schemes of Titarev and Toro [90 92 77J , the 
state vector was first expanded in a Taylor series in time at the element inter- 
face and second these time derivatives have been replaced by spatial deriva- 
tives using the so-called Cauchy-Kovalewski procedure that is based on a 
repeated use of the governing conservation law in differential form. The lead- 
ing state at the interface was computed by the exact Riemann solver applied 
to the boundary extrapolated left and right states. The higher order spatial 
derivatives at the interface have then been defined by solving linearized Rie- 
mann problems for the boundary extrapolated values of the derivatives, where 
the linearization was performed about the leading state. Though successful, 
the Cauchy-Kovalewski procedure, which is also equivalently called the Lax- 
Wendroff procedure [56] , turns out to be impracticable for complex systems of 
non-linear equations, and alternative methods need to be followed. Other ex- 
plicit one-step schemes in time that are based on the Lax-Wendroff procedure 
can be found, for example, in [721171 116 111431159] . 

While the Cauchy-Kovalewski procedure is based on the strong differential 
form of the PDE, an alternative has been proposed in [30.28J, where an 
element-local space-time Galerkin predictor is introduced that is based on 
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the weak integral form of the PDE. This approach is also able to account 
for stiff source terms. An overview of explicit one-step time discretization 
schemes can be found in [42J. The high order one-step schemes proposed in 
[301128] can be divided in three different steps. First, a high order WENO recon- 
struction |49ll91ll48ll8014QII51ll31l[32l! 104 93 197ll70j from the given cell averages ii 
is performed at time i n , but also other nonlinear high order reconstructions 
are possible, e.g. [46 8 4lfT1[2T] . Let the resulting piecewise polynomial solution 
be denoted by = w^(x, t n ) in the following. Second, the local space-time 
DG predictor [30.37J is applied for the time evolution of the reconstructed 
polynomials w^, the result of which are piecewise space^time polynomials de- 
noted by c\h = q/ l (x,t). This allows a high order accurate computation of 
the numerical fluxes. Finally, the cell averages are updated in time with a 
one-step scheme according to the integral form of the conservation law. The 
numerical method just described has been successfully applied to a variety 
of physical systems, including stiff advection-diffusion-reaction problems [47], 
compressible magnetohydrodynamics [28113] and magnetic reconnection in an 
astrophysical context [103] . 

Along with high order numerical schemes, another frontier of numerical re- 
search is represented by the implementation of efficient AMR algorithms, al- 
lowing for numerical simulations of very complex and highly dynamical struc- 
tures that require an adaptive refinement of the grids in specific flow regions. 
Starting from the pioneering AMR implementation by Berger et al. [T41I12II11] . 
who first introduced a patched-based block-structured AMR for finite differ- 
ence methods, several codes have been developed over the years providing 
AMR infrastructures in combination with different numerical techniques. We 
recall here, to list but a few, the AMRCLAW package [T5l[T3ll8] based on the 
second order accurate wave-propagation algorithm; AstroBEAR [241TL8] . where 
a collection of TVD and piecewise parabolic methods are applied for the spa- 
tial reconstruction, while a MUSCL-Hancock predictor-corrector scheme is 
adopted for the temporal evolution; RAMSES [89J, using a tree-based data 
structure with a second order Godunov method; NIRVANA [106] . with a sec- 
ond order, directionally unsplit central-upwind scheme of Godunov type [54] . 
Further AMR schemes can also be found in, e.g. [4 2 52 661167], 

In the last few years, moreover, significant progress has been obtained in com- 
bining high order numerical methods with AMR techniques. It is worth men- 
tioning the fourth order finite volume AMR scheme with Runge-Kutta time 
integrator by Colella et al. [23j; the PLUTO- AMR code by Mignone et al. [63], 
specifically devised for astrophysical applications, which includes a Corner- 
Transport-Upwind scheme for time updating and a third-order WENO spatial 
reconstruction; the spectral WENO schemes by Burger et al. [IT], combined 
with a third-order TVD Runge-Kutta method, meant for the study of the 
sedimentation of a polydisperse suspension. Finite difference WENO schemes 
together with AMR have been considered very recently in [79] . 
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Common in most of the above AMR methods is the use of the method of lines 
(MOL) together with a higher order (TVD) Runge-Kutta scheme in time 
[82ll83145j . In this paper we improve with respect to the approaches presented 
so far by providing the first tree-based 'cell-by-cell' AMR algorithm in combi- 
nation with a finite volume ADER-WENO one- step time update scheme. The 
use of high order one-step schemes in time allows directly for time-accurate 
local time stepping in a straight forward manner and has already been suc- 
cessfully applied in the context of high order Discontinuous Galerkin schemes 
with time accurate local time stepping [33118816 11143] . in multi-block domain 
decomposition methods [98] and finite volume and DG schemes based on static 
mesh adaptation [20] . 

The structure of the paper is the following. In Sect. [2] we provide a description 
of the ADER approach. Sect. [3] is specifically devoted to the description of 
the Adaptive Mesh Refinement infrastructure. In Sect. [| we present a variety 
of numerical tests for which two different systems of hyperbolic equations 
in conservative form have been considered, namely the classical Euler and 
magnetohydrodynamics equations. Finally, Sect. [5] contains the conclusions of 
our work. 



2 Numerical Method 

In this section, the numerical method is first presented for purely regular 
Cartesian meshes. Due to the one-step nature of our high order finite volume 
schemes, all basic ingredients can later be transferred also to AMR grids with 
only few modifications. 

2. 1 The finite volume scheme 

We consider hyperbolic systems of balance laws in Cartesian coordinates 

where u is the vector of conserved quantities, while f (u), g(u) and h(u) are the 
physical flux vectors along the x, y and z directions, respectively. As usual, 
we define the three-dimensional control volume on a regular Cartesian grid 
as I ijk = x [ yj _i;y j+ i] x [z k _i;z k +i], with Ax t = x l+ i - x t _i, 

Ayj = i —yj_i : Azk = z k+ i — z k _i, and, in addition, the spacetime control 

volume Iijk = Iijk x \k n -> t n + At] , where At = £ n+1 — t n . On adaptive meshes as 
used in this article, it is convenient to address each three-dimensional control 
volume with a unique mono-index m, like in the case of unstructured meshes, 
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with 1 < m < A^ceiis and ATceiis being the total number of cells at any given 
time (see Sect. 3.1). Therefore, in the rest of the paper we will also use the 
symbol C m to denote a control volume 1^. The numerical scheme, however, 
is more conveniently written using the indices z, j and k. The integration over 
lijk provides the standard finite volume discretization 

_ n+ i _ - n _ At / \ ^ At / ^ \ 

At 



Az k 

where 



-h ijtk _i)+AtS ijk , (2) 



x i y . , i 2, i 

* - T q Jt n "'"To 

111 - 

il n — 



= ~Kx' i ~Ky~ j ~Kzi : I I J^V^ndzdydx (3) 



3 " x ;-l y -,-i Z k -l 



"2 fc ~2 

is the spatial average of the solution in the element 1^ at time i n , while 



(4) 

t n + l^ i+ l^ + l 

gij+|,fc = At Ax / / J s(cih(x,y j+ i : z,t) : ql;(x,y j+ i : z : t))dzdxdt, 

(5) 

= At Ax- Ay I I J K^^^y^k+^^^t^^y^k+^^dydxdt, 



x. i y . i 



(6) 



and 
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t n+l ^+1 ^+1 



At Axi Ayj Az k 

t n x i y i z i 

2^22 



J J J J S(cLh(x,y,z,t))dzdydxdt . (7) 



are the space-time averaged numerical fluxes and sources, respectively. In the 
integrals above, the symbol denotes the local space-time DG predictor 
solution illustrated in section |2.3| In this article, we will use two alternative 
numerical fluxes, the classical Rusanov flux [74,94] or a new variant of the 
Osher-Solomon flux proposed in [36J. The Rusanov flux, also called local Lax- 
Friedrichs flux in literature, reads 

f = \ (f (<£) + f (<£)) " \\****\ (<tf - <fc) , (8) 
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where |<s max | denotes the maximum absolute value of the eigenvalues of the 
Jacobian matrix A = df/du. The Osher-type flux proposed in [36] reads 

f = \ (f(q^) + f (<tf)) " ^ (jf 1 \M^))\ds) (q£ - qu) , (9) 

with the straight-line segment path 

^(s) = q^ + s(q+-q^), < s < 1, (10) 

which connects the left and right state with each other in phase space. The 
path integral in Eqn. ^ is evaluated numerically using Gauss-Legendre quadra- 
ture with at least two quadrature points [36]. The numerical fluxes g and h 
are computed in the same way as the flux f . 

2.2 WENO reconstruction 

Since the method used here differs to some extent from the original WENO 
scheme of Jiang and Shu [19], some more details are given in the following. For 
each spatial dimension, we consider a nodal basis of polynomials of degree M 
rescaled on the unit interval I = [0; 1]. We remind that such a basis consists of 
M+l Lagrange interpolating polynomials of maximum degree M, {^/(A)}^ 1 , 
associated to the M + l Gauss-Legendre nodes {Xk}^ 1 in the interval [0; 1] 
and with the standard property that 

il>i(Xk) = Sik J,fc=l,2,...,M + l. (11) 

Here, 5ik is the usual Kronecker symbol. This choice produces by construction 
an orthogonal basis over [0; 1]. Moreover, it has the additional advantage that 
the data is immediately available in the Gaussian quadrature points anytime it 
is necessary to compute an integral over [0; 1]. The reconstruction is performed 
for each cell 1^ on a set of one-dimensional reconstruction stencils, which are 
given for each Cartesian direction by 

i+R j+R k+R 

Sijk = U left' ^ijk = U ^e/c, S*£ = IJ hjei (12) 
e=i—L e=j—L e=k—L 

where L = L(M, s) and R = i?(M, s) are the order and stencil dependent 
spatial extension of the stencil to the left and to the right, respectively. Odd 
order schemes (even polynomial degrees M) always adopt three stencils, one 
central stencil (s — 1, L — R — M/2), one fully left-sided stencil (s = 2, 
L = M, R = 0) and one fully right-sided stencil (s = 3, L = 0, R = M). Even 
order schemes (odd polynomial degree M) always adopt four stencils, two of 
which are central (s = 0, L =floor(M/2) + 1, R =floor(M/2)) and (s = 1, 
L =floor(M/2), R =floor(M/2) + 1), while the remaining two are again given 
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by the fully left-sided and by the fully right-sided stencil, respectively, as 
defined before. The total amount of cells of the stencil is the same as that of 
the order of the scheme, namely M + 1. After introducing a set of reference 
coordinates for each element 1^ given by 

x = x i _i+£Ax i , y = x-i+r]Ay j , z = x k _i + (Az k , (13) 

2 •* 2 2 

our reconstruction algorithm works according to a dimension by dimension 
fashion and can be described in the following steps: 



2.2.1 Reconstruction in x direction 

The reconstruction polynomial for each candidate stencil for element 1^ along 
the first coordinate direction x is written in terms of the basis function 
asH 

M 

w^ x (x, t n ) = MO™?jk, P Mt)™ijk, P • ( 14 ) 

p=0 

Integral conservation on all elements of the stencil then yields the following 
linear algebraic system for the unknown coefficients p of the reconstruction 
polynomial w^' x (x,t n ) in element 1^ 

^- f e+h *p{Z{x))*W<* dx = u n ejk , VI ejk e Stf k . (15) 

L\X e Jx i 

2 

For regular Cartesian meshes, the coefficients of the linear algebraic system 
above depend only on the choice of the basis functions, hence the system can 
be conveniently solved for the unknown coefficients by precomputing the 

inverse of the coefficient matrix, which can be done once and for all for each 
stencil on the reference element. Once the reconstruction has been performed 
for each of the stencils relative to the element 7^, we finally construct a data- 
dependent nonlinear combination of the polynomials obtained for each stencil, 
i.e. 

N s 

w x h (x : t n ) = ^(tfw^, with w^ p = £ u s w^ (16) 

s=l 

where the number of stencils is N s = 3 or N s = 4, depending on M being even 
or odd, respectively. The nonlinear weights are given by the usual relations 

m 

uj s = , oj s = — F , (17) 



1 Throughout this paper we use the Einstein summation convention, implying sum- 
mation over indices appearing twice, although there is no need to distinguish among 
covariant and contra-variant indices. 



7 



where the oscillation indicator a s is 



a s = EpmW^w^ , (18) 
and requires the computation of the oscillation indicator matrix [30j 



u 

In our implementation we have adopted A s = 1 for the one-sided stencils and 
A = 10 5 for the central stencils. Moreover, we use e = 10 -14 and r = 8. 

We stress that the resulting reconstruction polynomial w^(x,i n ) is only a 
polynomial in x direction, but still an average in the y and z direction, re- 
spectively. Hence, the reconstruction algorithm just described before can be 
applied again to the remaining two directions, as described below. 



2.2.2 Reconstruction in y direction 

When the reconstruction algorithm is applied along the second direction 



the steps from (TT4J) to (16) are repeated for each degree of freedom w^- fc . 
More precisely, wenave 



w^(x,y,n = MOMv)^f k , pq - (20) 
Integral conservation is now applied for each degree of freedom in x direction 
on all elements of the stencil S-jl in y direction, since the polynomial is still 
an average in the y direction. This yields 

XT / 6 3 Mv(y))™?j 8 k,pqd>y = w£ fc \/I iek G S-f k . (21) 

^Ve J y P _i 



' 2 



Again, the essentially non-oscillatory property of the reconstruction polyno- 
mial is assured using the nonlinear weighting of the individual reconstruction 
polynomials as 



N s 



W V h (x, y, t n ) = M&MV^ijkW With ™ijk*q = E ^ijkW ( 22 ) 



.s=l 



with uj s and a s defined as above in Eqn. (17) and (18). 



2.2.3 Reconstruction in z direction 

Finally, when the reconstruction algorithm is applied along the last direction 



z, the steps from (14) to (16) are repeated for the (M + 1) degrees of freedom 
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of the polynomial already reconstructed along x and y. We therefore have 



w£*(x, y, z, t n ) = ^ p (£)il> q (v)MC)™?jk,pqr ■ ( 23 ) 
with the integral conservation written as above, 

-!- f e+ * ^(CW)w^ dz = wZ k , pq , V/« e € fig. (24) 

•'z l 

e 2 

The final three-dimensional WENO polynomial is then given by 

w fc (x, t n ) = MOMv)MO^ jk , mr , (25) 

with 

N s 

™ijk,pqr = ^Z^s^ijk^pqr- ( 26 ) 
8=1 

We stress that the resulting polynomial WENO reconstruction produces entire 
polynomials, although they are given under the form of a nodal basis, but this 
is just a technical issue. Any other set of basis functions could have been 
chosen as well. This is the main difference with respect to the original optimal 
WENO scheme of Jiang & Shu [49] , which produces point values at the element 
interfaces. The drawback of our present method is of course that the method 
is not optimal, in the sense that to obtain a given order of accuracy the total 
stencil of the present method is larger than the one of the optimal WENO. 

2.3 Local space-time DG predictor 

Once a high order polynomial in space has been reconstructed for each cell, 
it is necessary to evolve it in time in order to compute the fluxes and sources 
according to Eqs. Q-([7]). Recall that the arguments of the numerical fluxes 
Q-([6]) are the yet undefined space-time polynomials q^. The strategy of a high 
order time evolution follows the spirit of the original MUSCL scheme of van 
Leer [100] as well as the one of the ENO scheme [46] . It can furthermore also be 
interpreted as an approximate solution of the generalized Riemann problem at 
the cell interfaces as done in the ADER approach [92,96], see fT9ll65] . However, 
the MUSCL scheme, the ENO method as well as the ADER approach obtain 
all the higher order in time through the use of Taylor series and a (repeated) 
use of the governing conservation law in strong differential form, which for 
higher than second order and for complicated PDE systems may quickly be- 
come very cumbersome or even unfeasible. Starting from the work by Dumbser 
et al. [30.28J an alternative time evolution has been proposed, which is able to 
deal with general nonlinear conservation laws, thus providing a more flexible 
one-step scheme compared to the ones based on Taylor expansions. The new 
method relies on a weak integral formulation of the governing PDE in space- 
time using an element-local space-time Galerkin predictor method. For this 
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time evolution strategy, no algebraic manipulations of space and time deriva- 
tives are necessary, but all that is required is a point-wise evaluation of fluxes 
and source terms. The result of the local space-time Galerkin predictor are 
the high order space-time polynomials that are needed for the evaluation 
of the numerical fluxes according to Eqs. @ ([7]). In the following we briefly 
illustrate the method, addressing to [30 28 371H7] for more details. 

In addition to the spatial reference coordinates £ = (£, 77, () already defined 



by Eqn. (13), we also introduce a reference time coordinate r as t = t n + rAt, 
spanning the unit interval [0, 1]. As a result, the governing PDE ([!]) can be 
rewritten in compact form as 

On df* <9g* <9h* 0# 

a7 + aTftT + lK =s (2T) 

with 

f* = ^f, g* = ^g, h* = ^h, S* = AtS. (28) 
Axi Ayj Az k 

We then introduce the space-time basis functions r), which are piecewise 
space-time polynomials of degree M, using the multi-index p = (p, g, r, s). The 
9p are given by a tensor-product of the basis functions ipi already used before 
in the reconstruction procedure, hence 

p (S,t) = MMMMOMt)- (29) 



Multiplication of the PDE ( |27| ) with the space-time test functions 9 q and 
integrating over the space-time reference control volume [0; l] 4 yields 



1 1 ifr + ^ + % + ^ ~ S " I * = ° (30) 



± ± ± ± 

////' 



Integration of the first term that contains the time derivative by parts yields 

111 liii/— \ 

a 



000 0000 
1111 



J J |# q (ei)u(ei)<wc- //// (^qj nd^dcdr 

d^drjd^dr 

JJJ q (CO)w h (tt n )dCdr ] dC. (31) 



+ 



0000 
1 1 1 



(df* dg* dh* N 







In the following, we denote the discrete space-time solution by q^, for which 
we make the following ansatz 

q/i = qfc(£, r) = 9 p (£, r) q p , (32) 
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with the yet unknown degrees of freedom of the space-time polynomial q p = 
Qpqrs We use the same representation for the fluxes and source terms, hence 



f,* = p f p *, gj = K = e p h;, 



q* n q* 



(33) 



Due to the nodal approach, the above degrees of freedom for the fluxes and 
source terms are simply the point-wise evaluation of the physical fluxes and 
source terms, hence 

t = r (q p ) , g! = g* (q p ) , h! = h* (q p ) , S: = S* (q p ) . (34) 



Inserting Eqns. (32) and (p3| into (31J) yields 
ill 



1111 



/ / /^ q (tl)^(tl)qp«- /// J (^eAe^d^drjdCdr 





1111 







+ 



(d % d d % 



III! 



111 

/ / Je^wn&nd^dc 



d^dr]d(dr 



(35) 







The above weak form is an element-local nonlinear algebraic equation system 
for the unknown coefficients q p . The initial condition is included in a weak 
sense by the integral on the right hand side, where the reconstructed solution 
w h(£,t n ) is given by Eqn. (25). After introducing the integrals 



ill 



1111 



K 



/ / fo q (&l)O p (t,l)d£- Jill (^oAOpdtdr, (36) 



ooo 



oooo 
1111 



K = K' K ? P ' K i) = //// e ^e p d^r, 



1111 

M qp = I J J I e,e p d£dT, 



and 





111 



(37) 
(38) 

(39) 



where d£ = dt;dr]d( one can rewrite the above system in compact matrix- 
vector form as 

K; p q p + K« p • I; + K;g p * + Kj p h; = M qp S p * + F; m <, (40) 

with the spatial multi-index m = (fc, /, m) and ^ m (£) = ipkiO^ii^^PmiC) - Due 
to the tensor-product nature of the basis functions and the control volumes, 



K = JJ /ft|(&0)Wa(€)#, 
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the above matrices are all very sparse block-matrices, where all sub-blocks 
contain purely one-dimensional integrals. Hence, the product of the matrices 
with the vectors of degrees of freedom can be efficiently implemented in a 



dimension-by-dimension manner. Equation (40) is conveniently solved by the 



following iterative scheme, introduced in [281137] . 

^qp^p 1VA qP^p — r qm W m ^qp T p ^qp&p ^qp^ V 4J J 

for which an efficient second-order MUSCL-type initial guess has been pro- 
posed in [47]. Moreover, the matrix Kj p can be easily inverted once and for all 
on the reference element by inverting the one-dimensional sub-blocks in time 
that are associated with each Gaussian quadrature point in space. 



3 Adaptive Mesh Refinement 



There are two major strategies for implementing an AMR algorithm. The first 
one adopts nested arrays of logically rectangular grid patches, according to the 
original Berger-Colella-Oliger approach [T4lll2llllj . The second strategy is re- 
ferred to as the 'cell- by-cell' refinement, the implementation of which is some- 
how similar to the one of finite volume schemes on unstructured meshes [53] . 
In our work we have followed the second choice, since it can be easily imple- 
mented using a tree-type data structure and is slightly more general than the 
first one. 

As already mentioned at the beginning of the previous section, the three in- 
gredients of our high order one-step finite volume schemes (reconstruction, 
time evolution, finite volume update) can be transferred in a rather straight- 
forward manner to space-time adaptive grids. In particular, the element-local 
space-time DG predictor on regular Cartesian grids is identical to the one 
on AMR grids, since for the predictor there is no need to exchange informa- 
tion with neighbor elements. Hence, even if two adjacent cells are on different 
levels of grid refinement this will not alter the local space-time DG predic- 
tor scheme at all. The reconstruction procedure obviously needs cell averages 
from neighbor elements. Therefore, in our AMR implementation the case of 
adjacent cells with different levels of refinement is handled in such a way that 
a real cell of grid level £ is always surrounded by a sufficiently thick layer of 
virtual (ghost) cells on the same layer. The cell averages of the virtual cells 
are obtained from the real cells on the same location either by averaging or 
projection, depending whether the virtual cell is on a coarser or finer level 
of refinement. Since for high order WENO schemes the total effective stencil 
needed for obtaining the reconstruction polynomial within a spatial control 
volume Iijk grows with the degree of the reconstruction polynomial M, the 
layer of ghost cells must always have at least a thickness of M ghost cells. In 
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order to keep the information needed for reconstruction local on the coarser 
grid level, the grid refinement factor t adopted in the scheme [see definition 
(42) below] must satisfy t > M. The fact that active cells are surrounded by 
a layer of ghost cells can of course also be interpreted as use of micro-patches. 
The virtual ghost cells are logically connected neighbors of the real cells in or- 
der to allow the uniform Cartesian reconstruction procedure to be performed 
exactly as described in the previous section. The only real difference between 
uniform Cartesian grid and AMR grid can be found in the computation of 
the fluxes at edges with two adjacent cells of different level of refinement. To 
keep the algorithmic complexity as low as possible, in our approach the level 
of refinement of two adjacent grid cells may differ by at most one. 



3.1 AMR implementation 



We have developed a cell-by-cell AMR technique in which the computational 
domain is discretized with a uniform Cartesian grid at the coarsest level. We 
use Cq to denote this initial grid on the coarsest level of refinement (t = 0), 
while Ct indicates the union of all elements up to level £. Already at time t — 0, 



the refinement criterion (see Sect. 3.1.1) is applied to the initial condition, 



thus producing a hierarchy of refinement levels up to a prescribed maximum 
level of refinement ^ max , i.e. < ^ < ^ max - In the rest of this section we will 
use the following terminology: With children we intend the cells on the next 
refinement level £ + 1 contained in a cell C m of level £ after its refinement. 
For the children cells the original cell C m is denoted as their mother cell. The 
Neumann neighbors M m of a cell C m are the neighbor cells that share a common 
face with cell C m . There are 2d Neumann neighbors in d space dimensions. 
For example, the numerical fluxes (|4])-([6]) are evaluated between Neumann 
neighbors. The Voronoi neighbors V m of a cell C m are cells which share common 
nodes, hence the Neumann neighbors are a subset of the Voronoi neighbors. 
There are 3 d — 1 Voronoi neighbors in d space dimensions. 

All over the simulation, the following general rules have to be met: 

(1) Whenever a cell of the level £ is refined, it is subdivided into an integer 
number t of finer cells along each direction, such that 

Ax £ = xAx w Ay i = xAy £+1 Az £ = xAz £+u (42) 

and also the time steps are chosen locally on each level so that 

At £ = xAt £+1 . (43) 

As a result, each mother cell generates x d children cells in d space dimen- 
sions. Since the high order WENO reconstruction needs information from 
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more cells than just the direct neighbors, the following condition must 
always hold to keep reconstruction local on the coarser grid level: 

t> M. (44) 

(2) Each cell C m , at any level of refinement, has one among three possible 
status flags. The status flag is denoted by a in the following. A cell is 
either an active cell (a = 0), and therefore has to be updated through the 
finite volume scheme; or, it is a virtual child cell (with status a = 1) and is 
updated by projection of the mother's high order space-time polynomial 
q^; or, finally, a virtual mother cell (a = — 1), updated by recursively 
averaging over all children from higher refinement levels. A virtual child 
cell has always an active mother cell, while a virtual mother cell has 
children with status a < 0. A virtual child cell cannot be further refined, 
unless it is first activated. 

(3) Any cell C m , at any level of refinement and with any status, is identified 
with a unique posit iv^] integer number m, with 1 < m < Nc e \\ s and 
A^ceiis being the total number of cells at any given time. ATceiis is of course 
a time-dependent quantity, increasing for any refinement operation, and 
decreasing for any recoarsening operation. 

(4) A refined cell with children of status a < cannot have any Voronoi 
neighbor without children. Therefore, as soon as a cell C m is refined and 
the generated children have status a = 0, or as soon as the existing virtual 
children of cell C m are activated, all Voronoi neighbors without children 
are virtually refined, i.e. they generate virtual children with status a = 1. 
The status of the regularly refined cell C m changes to a — — 1. 

(5) The levels of refinement of two cells that are Voronoi neighbors of each 
other can only differ by at most unity. Violation of this rule is avoided 
through appropriate activation of virtual children. Such a situation is 
schematically depicted in Fig. [T} 

(6) The maximum level of refinement is a prescribed level ^ max . Cells are not 
allowed to be refined beyond this level. 

(7) When a cell C m is recoarsened, its children cells are destroyed. However, 
this is not allowed if the cell C m contains children that contain themselves 
children of any status. Furthermore, if the cell C m is to be recoarsened 
and has a neighbor cell with active children, then the cell C m can only 
deactivate its children cells, changing them to virtual, but cannot destroy 
them. 



For each cell we store the indices of the Voronoi neighbors, as well as pointers 
to the mother and the first child of a cell. The indices of all x d children are 
obtained by the convention that all children of a cell have consecutive numbers. 
We furthermore store the status a of all elements. 



See Sect. 3.2 for the possibility of negative integer numbers assigned to those 



cells belonging to the ghost zone at the MPI border between two processors. 
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Fig. 1. Refinement mechanism involving two levels of refinement (£ < 2) in a 2D 
geometry. Solid lines denote active cells, while dashed lines denote virtual cells. Left 
panel: when the central cell Cj at level £ = is refined (solid-dark-red line), all its 
Voronoi neighbors Vj are virtually refined (dashed-light-red lines). Right panel: if the 
cell Ck at level £ = 1 on the left corner of Cj requires further refinement to level £ — 2 
(solid-dark-blue lines), the virtual children of the lower, left and lower left neighbors 
of Cj must be activated (to avoid violation of rule (5)), hence all their Voronoi 
neighbors without children must be virtually refined to level £ — 1. Furthermore, 
the neighbors of Ck must be virtually refined to level £ — 2 (dashed-light-blue lines). 

A simple ID example of the tree-structure used in our algorithm together 
with an illustration of the use of the status flag a is depicted in Fig. [2j 



3.1.1 Refinement criterion 

As a refinement criterion, we have adopted the same strategy described in [60] 
and later adopted by [41J and [63]. Such a criterion has two very desirable 
properties. The first one is that it is entirely local, thus avoiding any global 
computation. The second one is that it is based on the calculation of a second 
derivative error, thus avoiding unnecessary refinement in smooth regions, such 
as the wake of a rarefaction wave. In practice, a cell C m is marked for regular 
refinement (generation of children with status a = 0) if Xm > Xref ? while it is 
marked for recoarsening if Xm < Xrec, where 



Y.U&*ld*kdxif 

Zk,il(\d<S>/dx k \ i+1 + IdQ/dx^/Axt + e\{d*/dx k d Xl )\\n 2 ' 1 j 

The summation J2ki i s taken over the number of space dimension of the prob- 
lem in order to include the cross term derivatives. We stress that $ = $(u) 
can be any suitable indicator function of the conservative variables u. We have 
observed that the threshold values Xref and Xrec can be slightly model depen- 
dent. In most of our tests we have chosen Xref i n the range ^ [0.2, 0.25] and Xrec 
in the range ^ [0.05, 0.15]. Finally, the parameter s acts as a filter preventing 
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Fig. 2. One-dimensional example for the tree structure used in the present cell- by- 
cell implementation of AMR with r = 3 and t < 2. 

refinement in regions of small ripples and is given the value e = 0.01. 



3.1.2 Time accurate local timestepping 

Each of the refinement levels is advanced in time with its own local time- 
step Ati = tA^ + i, as already defined above. The use of time steps that are 
integer multiples of each other for each level is very convenient, but not strictly 
necessary, see e.g. the local time stepping schemes presented in [33ll88ll61j . 
where a different local time step is allowed for each element. Let us denote by 
t™ and the current and future times of the level L Then, the level scheduled 



for updating is the largest value of t that satisfies the update criterion [33j 

t n , +1 <t n +l 0<£<Cax, (46) 



where we define t 7 ^ 1 := £q +1 for convenience, so that also the scheduling of 



level I — is included in Eqn. (46). 



In other words, starting from the common initial time t = 0, the finest level 
of refinement ^ max is evolved first and performs a number of t sub-timesteps 
before the next coarser level ^ max — 1 performs its first time update. This 
procedure is then applied recursively and it implies a total amount of x £ sub- 
timesteps on each level to be performed in order to reach the time t$ +1 of the 
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coarsest level. 



The computation of numerical fluxes between two adjacent cells on different 
levels of refinement is rather straightforward thanks to the use of the local 
space-time predictor, which computes the predictor solution for each ele- 
ment after reconstruction and which is valid from time t™ to time In this 



way, and with the update criterion (46), a high order space-time polynomial 



is always available on both sides of an element interface for flux computation. 

To illustrate the procedure, we use the one-dimensional example depicted in 
Fig. [3j which corresponds to the case already shown in Fig. [2] before. All ele- 
ments start from a common time level an d the space-time predictor solution 
q/j has been computed in all elements. Then, the first level that satisfies the 



update criterion (46) is t = 2. Computing the fluxes for cell C 19 is exactly as 
on a uniform Cartesian mesh, since C\g has two active neighbors on the same 
grid level. The situation is different for C\% because its left real (active) neigh- 
bor is cell Cii, which is on a coarser grid level. Nevertheless, the numerical 
flux between C 18 and C\\ can be computed since the predictor solution is 
available for the necessary time interval [£ 2 ;£ 2 + Ai 2 ] in both cells. To make 
the approach conservative, the computed numerical flux between C n and C±$ 
is used to update Ci 8 , but it is at the same time also stored in a memory vari- 
able of the virtual neighbor cell C23 and will be used later to update the real 
element C\\. For convenience, the virtual neighbor cell C 23 can also be used to 
technically handle the hanging node in time (and in space for d > 1) by scal- 
ing the space-time polynomial q^ of the coarse element C\\ down to the finer 
grid level £ = 2. A similar situation occurs for element C 2 q and its right active 
neighbor C 7 . After the first time step of level I = 2, the elements C 18 — C 20 
are at time t$ + At 2 , while all other elements are still at time To perform 
reconstruction for cells C 18 — C 20 , the predictor solution q^ is projected from 
elements C\\ and C7 into the cell averages u of their virtual children at time 
£0 + At 2 . Now, reconstruction for cells Ci 8 — C 2 o can be performed exactly as for 
the uniform Cartesian case and subsequently the local space-time predictor 
can be carried out with the new initial data w^. For the time update in the 
next time interval [tg + Ai 2 ; t$ + 2Ai 2 ] the predictor solution in elements C\\ 
and C7 is still valid, hence the numerical flux can again be directly computed 
between cells as described before. This procedure is carried out until all cells of 
level £ = 2 reach the time t$ + tAt 2 . Then, their future time is ifi + (t+ l)Ai 2 , 



hence the update criterion (46) is no longer fulfilled for level £ = 2 and the 
next coarser level I = 1 can be updated. Cell C10 is only surrounded by active 
neighbors on the same level, hence the standard finite volume scheme on uni- 
form Cartesian mesh can be used. For cells C 9 and C 8 the previously described 
procedure of flux computation between fine and coarse cell applies. Cell C\\ 
now illustrates the last special case, namely the computation of the numerical 
flux between the coarse cell C\\ and fine cell C\%. Actually, the solution is par- 
ticularly simple. Due to the requirement of conservation, the coarse cell C\\ 
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must not compute any new flux at the interface with element Cis, since all the 
necessary fluxes have already been computed before on the finer level. It is 
just sufficient for cell C\\ to sum up the fluxes stored in the memory variable 
of its virtual child C23. In this way, the sum of time integrals of the fluxes on 
the left border of element C\% is equal to the time integral of the fluxes on the 
right border of C\\ in the time interval [£q; £q + Ati], which makes the method 
conservative. In the most general 3D case, the numerical flux computed at an 
interface containing elements of level t and £ + 1 reads 

f -+^ = t tt IJJ (47) 



with the integration intervals above defined as 



Ta = [*? + (ii - l)At m ; t n t + iiAf m ], 
^•j = + (jj - l)A^+i; + jjA^+i], 

^fcfc = + (fefe - l)A^ + i; + fcfcA^+i]. (48) 



The flux on the left hand side of (47) corresponds to the final averaged flux 
for the coarse grid cell, while the integrals on the right hand side of (47) are 
the integrated fluxes for each time step of each fine grid cell adjacent to the 
coarse cell. The practical implementation of Eqn. (47) is conveniently achieved 
by the sum over the memory variables, as described before. 



This completes the description of all cases to be treated by the AMR algorithm 
concerning flux computation across elements on different refinement levels. 



3.1.3 Projection 

Projection is the typical AMR operation, sometimes called " coarse-to-fine pro- 
longation", by which values an active mother assigns values to the virtual 
children (a = 1) at intermediate times via standard L 2 projection. For this 
purpose, the space-time polynomials can be conveniently evaluated at any 
time. This operation is needed for performing the reconstruction on the finer 
grid level at intermediate times. The projection operator for a cell C m on level 
t is simply given by evaluating the space-time polynomial of its mother at 
any given time as follows: 
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Fig. 3. One-dimensional example for the local timestepping with t = 3 and t < 2. 
The grid and the element numbers correspond to the ones depicted in Fig. [2j 

3.1.4 Averaging 

Averaging is another typical AMR operation by which a virtual mother cell 
(a = — 1) obtains its cell average by averaging recursively over the cell averages 
of all its children and their children at higher refinement levels. Let us denote 
the set of children of a cell C m by B m , then the averaging operator is given by 

™rn = ^ ^ ( 50 ) 

3.2 Overall efficiency and MPI parallelization of higher order AMR 

In the following we study quantitatively the overhead introduced by the high 
order one-step ADER-WENO AMR method proposed in this article. For this 
purpose, we report the CPU times needed for the simulation of the two- 
dimensional explosion problem discussed in more detail in section [4] for uni- 
form and AMR grids for second to fourth order ADER-WENO schemes. The 
detailed CPU time results are summarized for all cases in Tabled] and are nor- 
malized with respect to the standard second order scheme on uniform mesh. 
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The data refer to the average CPU time needed for the update of one single 
real element, which has been computed by dividing the total wallclock time 
needed for the simulation by the number of time updates of the active ele- 
ments (a = 0) contained in the domain. Hence, the results reported in Table 
[T] include the entire overhead necessary for the update, averaging and projec- 
tion of the virtual ghost cells needed in the AMR approach. In the table, we 
also report separately the total overhead introduced by the AMR approach 
in percent for convenience. The CPU times have been obtained on one single 
core of an Intel i7-2600 CPU with 3.4 GHz clock speed and 12 GB RAM. 

We furthermore have parallelized the three dimensional ADER-WENO code 
through the standard Message Passing Interface (MPI). Any AMR implemen- 
tation poses additional challenging problems to the parallelization task, which 
become manifest when refinement of a cell C m occurs at the MPI border be- 
tween two processors (see left panel of Fig. In this case, in fact, proper 
communication among the processors must be established in order to spread 
the knowledge about which cells must be either virtually refined or activated. 
For this purpose, each processor stores in its memory also MPI ghost-cells 
that are a copy of the true cells, managed by the adjacent processor. In the 
practical implementation, we have found convenient to assign a negative inte- 
ger number to each cell in the MPI ghost-zone, thus making the distinction 
with respect to real cells very transparent. When a cell in the domain of the 
processor CPUO, at the border with the domain of the processor CPU1, is re- 
fined (see right panel of Fig. [4]), CPUO informs CPU1 that (i) a number of real 
cells belonging to CPU1 must be virtually refined, and (ii) that one cell (the 
one at the border) must be virtually refined in the MPI-ghost zone of CPU1. 
This information is used by CPU1 to (virtually) refine its corresponding cells 
in its MPI-ghost zone. Before doing that, CPU1 must also check whether 
such cells have already received an instruction of virtual refinement internal 
to CPU1. The link between the true cells of CPUO and those belonging to 
the MPI-ghost zone of CPU1 is obtained via so-called exchange lists. The 
exchange lists are used for the MPI communication that is necessary during 
the adaptive mesh refinement procedure, as well as to exchange the informa- 
tion about the cell averages u and the space-time polynomials between the 
processors. We stress that our present MPI- AMR implementation does not yet 
provide dynamic load-balancing among processors, which is a rather complex 
topic that will be considered in the future. 

Table 1 

Assessment of the overall efficiency of high order one-step ADER-WENO schemes 
on space-time adaptive AMR grids (t = 4,^ = 2). Normalized average CPU time 
per real element update with respect to the second order scheme on uniform grid. 



Scheme order 


Uniform grid 


AMR grid 


Total AMR overhead 


02 


1.00 


1.15 


15 % 


03 


3.18 


3.82 


20 % 


OA 


8.64 


10.82 


25 % 



20 



CPU1 



CPUO 



+4- 



+4- 



ghostlcells I — L _ 
I L J_ 



J_ 1 . 



j — i-i — i 

J_ _L _L J 



n 
a 



cells 


-\-'r 


— j - j — 


— i - 1 — i 
— 1 - 1 — ' 






i i 








i i _ 
- i - 1- - 





















c 
o 



Fig. 4. Cell refinement at the MPI-border between two processors. Left panel: pro- 
cessors CPUO and CPU1 must exchange information about which cells are refined 
(solid red line) or virtually refined (dashed black line). Right panel: Each processor 
has a MPI-ghost zone of cells that are a copy of the true cells managed by the 
adjacent processor. 

4 Numerical Tests 



In all the following numerical test problems we have used the density as indi- 



cator function in (45), hence <fr(u) = p. 



4-1 Euler equations 



The first session of tests considers a sequence of applications for which the 
classical Euler equations of compressible gas dynamics are solved. In three 
space dimensions the vectors of the conserved variables u and of the fluxes f , 
g and h are given respectively by 
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pv 2 z +p 










[vy(E+ P ) ) 




K v z (E + p) j 



(51) 

where v X) v y and v z are the velocity components, p is the pressure, p is the 
mass density, E = p/{^ — 1) + p{v 2 x + v 2 y + v 2 z )/2 is the total energy density, 
while 7 is the adiabatic index. 



2D isentropic vortex. The first test considered is a two-dimensional con- 
verted isentropic vortex, see e.g. [48]. The computational domain is Q = 
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Table 2 

Numerical convergence results for the isentropic vortex test using the third and 
fourth order version of the one-step ADER-WENO finite volume scheme presented 
in this article. The error norms refer to the variable p (density) at the final time 
tf = 10. The asterisk * refers to a uniform grid. 
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[0; 10] x [0; 10] and the initial conditions are given by a perturbation added to 
a uniform mean flow 



(p,v x ,v y ,v z ,p) = {1 + Sp, 1 + Sv x ,l + Sv y , 0,1 + 5p), 



(52) 
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The perturbation ST in the temperature is 

6 2 ( 7 -l 



ST 



877T 2 



exp (1 — r ) , 



(53) 



(54) 



with r 2 = (x — 5) 2 + (y — 5) 2 , vortex strength e = 5 and adiabatic index 
7 = 1.4. The refinement factor adopted is t = 3. In Table [2] we have reported 
the results of the convergence tests, where we have used the third and fourth 
order version of the method. The convergence rates have been computed with 
respect to an initially uniform mesh, as proposed by Berger and Oliger in [14]. 
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Fig. 5. Isentropic- Vortex test at the final time t = 10. Left panel: Contour plot of 
the mass density. Right panel: zoom into the AMR grid. Two levels of refinement 
have been adopted (^ max = 2). 



Interacting blast waves. This test, originally proposed by |l()2j . has by 
now become a classical problem in computational fluid dynamics and consists 
in the interaction of blast waves with initial conditions given by 



(1.0, 0.0, 10 3 ) if -0.5 < x < -0.4, 
(1.0,0.0, 10- 2 ) if -0.4 < x < 0.4, 
(1.0,0.0, 10 2 ) if 0.4 < x < 0.5. 



(55) 



Although one dimensional, we have evolved this problem in two spatial di- 
mensions over the domain [—0.5,0.5] x [—0.5,0.5], using reflecting boundary 
conditions in x direction and periodic boundary conditions along the y di- 
rection. The adiabatic index has been chosen as 7 = 1.4. Fig. [6] shows the 
results of our test, where we have used two levels of refinement from an origi- 
nal uniform grid with 300 x 10 cells, and adopting a third order ADER-WENO 
scheme. The left panel shows the solution at the time when the two waves hit 
each other from opposite directions, producing a very strong density peak. The 
right panel, on the other hand, shows the solution at the final time after the 
waves have crossed each other. A reference solution is also reported, obtained 
with a traditional finite difference TVD method using 3600 grid-points. 



Explosion problems in two and three space dimensions. In this 
test, proposed in [92 ,94] , we solve the Euler equations on the computational 
domain Q = [—1; l] d , where d denotes the number of space dimensions. The 
initial flow variables take constant values for r < R and for r > i?, separated 
by a cylindrical or spherical discontinuity, respectively. Therefore the initial 



23 




Fig. 6. Interacting Blast Waves Test. Profile of the mass density at time t — 0.028 
(left panel) and at the final time t — 0.038 (right panel). Two levels of refinement 
have been adopted from an initial grid 300 x 10. 

condition is given by 

{ix; if r < R, 
(56) 
u if r > R. 

Here, R = 0.4 denotes the radius of the initial discontinuity, x is the vector 
of spatial coordinates with the radial coordinate r = Vr. and u are the 
inner and outer states, respectively, listed in detail in Table [3j The adiabatic 
index of the ideal-gas equation of state has been set to 7 = 1.4. Due to the 
symmetry of the problem, which is cylindrical in the two-dimensional case and 
spherical in the three-dimensional case, the solution can be compared with an 
equivalent one dimensional problem in radial direction r, see [94J: 

Table 3 

Inner and outer initial states for the multidimensional explosion test problems. The 
last column reports the final simulation time t e . 
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(57) 



where u is the radial velocity component. The ID reference solution has been 
computed by a classical second order TVD finite volume scheme on a very fine 
mesh composed of 10000 grid zones and using the Osher-type flux proposed 
in [36]. The two-dimensional AMR simulations have been carried out with a 
fourth order ADER-WENO scheme on a level zero grid with 34 x 34 control 
volumes, using the Osher flux ^ and t = 3 and ^ max = 2. This leads to 
an equivalent resolution on a uniform fine grid of 306 x 306 = 93636 points. 
Fig. [7] shows a 3D plot of the density distribution obtained for the cylindrical 
explosion case, as well as the AMR grid configuration at the final time t = 0.25. 
Fig. [8] shows the results obtained on a one-dimensional cut along the x-axis, 
together with the ID reference solution according to (57). The cut is performed 



on equidistant points, evaluating locally the reconstruction polynomials w^. 
For comparison, we also show the results obtained on the uniform fine grid 
of 306 x 306 grid points that corresponds to the finest AMR grid level. Both 
simulations agree very well with the ID reference solution. Furthermore one 
can note only very little differences in the numerical results obtained with 
AMR and without AMR, i.e. on the uniform grid. However, the simulation on 
the AMR grid took only 4 minutes on 4 cores of an Intel i7-2600 CPU with 3.4 
GHz clock speed and 12 GB of RAM, while the fine uniform grid computation 
needed 14 minutes on 4 cores on the same machine, hence it took 3.5 times 
longer. This clearly confirms that even though the use of AMR adds a certain 
overhead of about 25 % to the fourth order finite volume scheme, according 
to Table [TJ the use of space-time adaptive meshes can significantly speed 
up multidimensional computations also for higher order schemes. The total 
number of AMR grid cells present at the final time was 28036, compared to 
the 93636 cells of the uniform grid. 



Similarly, Fig. 10 illustrates the solution for the three-dimensional problem, 
for which again a level zero mesh with 34 x 34 x 34 cells has been used, 
together with t = 3 and ^ max = 2. This corresponds to an equivalent fine 
grid resolution of 306 3 = 28, 652, 616 cells. In the three-dimensional case, a 
third order ADER-WENO scheme has been employed based on the Osher 
flux Q. The final AMR grid at time t = 0.25 contains 9,079,984 elements 
and is depicted in Fig. [9| The simulation took 7.5 hours on 8 CPU cores of an 
AMD Opteron 6272 cluster with 2.1 GHz clock speed and 256 GB of RAM. 
One-dimensional cuts (along the positive x-axis) through the reconstruction 



polynomials are shown on equidistant points in Fig. 10 together with the 



ID reference solution according to (57). An excellent agreement is observed. 
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Fig. 7. Explosion test in two space dimensions. Density distribution at time t = 0.25 
obtained with a fourth order ADER-WENO scheme (left) and corresponding final 
AMR grid configuration (right). 



1 D reference sol 
ADER-WENO OA 
ADER-WENO OA 



1 D reference sol 
ADER-WENO OA 
ADER-WENO OA 




Fig. 8. Explosion test in two space dimensions. One-dimensional cut along the pos- 
itive x-axis through the fourth order ADER-WENO solution obtained on the AMR 
grid for density (left) and velocity (right). The solution computed on a uniform fine 
mesh corresponding to the finest AMR grid level is also shown. 

Double Mach reflection problem. A classical test problem that con- 
tains simultaneously strong shock waves, contact waves and shear waves is 
represented by the double Mach reflection test [102], whose initial conditions 
are given by a right-moving shock wave at shock Mach number M = 10, 
intersecting the x— axis at x = 1/6 with an inclination angle of a = 60°. 
The computational domain is Q = [0;3] x [0;1]. On the left, top and right 
boundary, the exact solution is prescribed, while a reflective wall boundary 
is imposed on the bottom. This test problem is frequently used in the lit- 
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Fig. 9. Explosion test in three space dimensions. AMR grid structure at time t = 0.25 
and density contour colors. 



erature on high order WENO and Discontinuous Galerkin schemes, see e.g. 
[49114 6 6 78 22 8 TU921I71 II 1 051132] , where many reference solutions for this prob- 
lem can be found. 

The initial condition for this problem is given by the Rankine-Hugoniot con- 
ditions as follows: 



(p,v x ,v y ,v z ,p) 0,?/, 0) = < 



(8.0, 8.25 cos(a), 8.25 sin(a),0., 116.5) if x f < 0.0, 

(1.4, 0.0, 0.0., 0., 1.0) if x' > 0.0, 

(58) 

with x' = (x — 1/6) cos(a) — ys'm(a). 

The ratio of specific heats is chosen as 7 = 1.4. The problem is solved with 
a third order ADER-WENO scheme using the Rusanov flux ([8]). We use a 
mesh on the coarsest level consisting of only 150 x 50 elements, together with 
a refine factor of t = 4 and a maximum refinement level of ^ max = 2. On the 
finest level, this corresponds to an effective resolution of 2400 x 800 control 
volumes. 



The results for the density (31 equidistant contour levels from 1.5 to 22.5) are 
depicted at time t = 0.2 in Figure [TT] together with the final AMR grid. A 



zoom of the solution and the mesh is shown in Fig. 12 
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Fig. 10. Explosion test in three space dimensions. From top left to bottom right: 
comparison of the ID reference solution with the numerical solution obtained with a 
third order ADER-WENO scheme on AMR mesh at time t = 0.25. One-dimensional 
cuts along the x-axis are shown for density, velocity, pressure and internal energy. 

The shear waves present in this test are subject to the classical Kelvin- 
Helmholtz instability and therefore tend to roll up. Since there is no physical 
viscosity in the compressible Euler equations solved here, the developed small- 
scale flow features are purely governed by numerical viscosity. However, the 
amount of roll-up is a good qualitative indicator of the amount of numerical 
viscosity since more roll up indicates less numerical viscosity introduced by 
the scheme. The results obtained with the present ADER-WENO scheme on 
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X 



Fig. 11. Double Mach reflection problem at time t = 0.2. Top: equidistant density 
contour lines (contour spacing Ap = 0.5). Bottom: AMR grid with two levels of grid 
refinement. 




Fig. 12. Zoom into the double Mach reflection problem at time t = 0.2. Left: 
equidistant density contour lines (contour spacing Ap = 0.5). Right: AMR grid 
with two levels of grid refinement. 

space-time adaptive Cartesian meshes is in good qualitative agreement with 
other published results for this test problem. 
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Forward facing step. Another classical test problem for high resolution 
shock-capturing finite volume scheme consists in the forward facing step prob- 
lem, also called the Mach 3 wind tunnel test. It has also been proposed origi- 
nally in [102J. The computational domain is given by Q = [0; 3] x [0; 1]\[0.6; 3] x 
[0; 0.2] and the initial condition is a uniform flow at Mach number M = 3 
moving to the right. In particular, we use p(x,y, 0) = 1, p(x,?/,0) = I/7, 
v x {x ) y, 0) = 3 and v y — v z = 0. The ratio of specific heats is set to 7 = 1.4. 
Simulations are carried out until t = 2.5. Reflective boundary conditions are 
applied on the upper and lower boundary of the domain and inflow/outflow 
boundary conditions are applied at the entrance/exit. At the corner of the step, 
there is a singularity, which is properly resolved with the third order ADER- 
WENO scheme using adaptive mesh refinement. The mesh on the coarsest 
level contains 150 x 50 control volumes. We use t = 4 and ^ max = 2, hence on 
the finest level this corresponds to an equivalent resolution of 2400 x 800. The 
computational results obtained with the third order ADER-WENO method as 
well as a sketch of the final AMR mesh are depicted in Fig. [13J For comparison, 
also a second order simulation is shown. One can clearly observe that the third 
order scheme provides a much better resolution of the physical instability and 
roll up of the contact line compared to the standard second order scheme. This 
indicates that even in the context of space-time adaptive mesh refinement, the 
use of higher order schemes may be appropriate to enhance resolution and to 
reduce numerical viscosity for small scale turbulent structures. 



2D Riemann problems. A large set of two-dimensional Riemann prob- 
lems has been cataloged in [55j. The computational domain is Q = [— 0.5; 0.5] x 
[—0.5; 0.5] and the initial conditions are given by 



u(x,2/,0) = < 



Ui 


if 


X 


> Ay > 0, 


u 2 


if 


X 


< Ay > 0, 


u 3 


if 


X 


<0At/<0, 


u 4 


if 


X 


> Ay < 0. 



(59) 



The initial conditions and the final simulation time tf for the four configu- 
rations presented in this article are listed in Table [4} In all cases 7 = 1.4. 
The simulations are carried out with a third order one-step ADER WENO 
scheme using a level zero grid of 50 x 50 elements. The computational results 
together with the final AMR grids are depicted in Fig. [T4j We can note a good 
agreement with the reference solution published in [55j . 



A co-rotating vortex pair. This multi-scale problem from aeroacoustics 
solved here differs from the previous ones in several aspects. First, it is a low 
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Fig. 13. Forward facing step problem. AMR grid (top) and density contours (center) 
obtained with the third order ADER-WENO scheme at time t = 2.5. One can clearly 
observe the roll up of the slip lines. For comparison, also a second order solution is 
shown (bottom). 
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Fig. 14. Two-dimensional Riemann problems solved with third order ADER-WENO 
schemes. Density contour lines (left) and AMR grid at the final time (right). 
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Mach number problem without shock waves. Second, it is of true multi-scale 
nature. The problem is an academic example of sound generation mecha- 
nisms in aeroacoustics and is taken from |64157II25H75II98II73| . It consists of two 
isentropic vortices with characteristic size r c (vortex core radius) that rotate 
around each other, thus generating sound waves with a wavelength that is 
about two orders of magnitude larger than the size of the vortices themselves. 
While the sound waves are very smooth, the vortices contain strong gradients 
in velocity and pressure as one approaches the vortex core. The accurate prop- 
agation of sound waves with little amplitude and phase errors in large domains 
is the major challenge in computational aeroacoustics. For such problems, typ- 
ically special high order schemes are employed, e.g. [58. 85, 86, 76, 34], since the 
conventional high resolution TVD schemes that are typically used in classi- 
cal AMR codes are too diffusive and lose accuracy at local extrema. Since 
such extrema regularly occur in acoustic wave propagation problems the use 
of high order WENO schemes, as the ones used in this paper, that are able to 
handle strong gradients without degenerating at local extrema, seems to be 
appropriate. 

The initial condition for the velocity field is given by the superposition of the 
velocity fields induced by two potential vortices. The complex potential w of 
the rotating vortex pair is given by 

W ( z ,t) = ^-\n(z 2 -b 2 ), (60) 

Table 4 

Initial conditions for the two-dimensional Riemann problems. 



Problem RP1 (Configuration 3 in [55J), t f = 0.25 





x<0 


x > 




p u v p 


p u v p 


v >o 
y <o 


0.5323 1.206 0.0 0.3 
0.138 1.206 1.206 0.029 


1.5 0.0 0.0 1.5 
0.5323 0.0 1.206 0.3 


Problem RP2 (Configuration 4 in [55J), t f = 0.25 




x < 


x > 




p u v p 


p u v p 


v >o 
y <o 


0.5065 0.8939 0.0 0.35 
1.1 0.8939 0.8939 1.1 


1.1 0.0 0.0 1.1 

0.5065 0.0 0.8939 0.35 


Problem RP3 (Configuration 6 in [55]), tf — 0.30 




x < 


x > 




p u v p 


p u v p 


y >o 
y <o 


2.0 0.75 0.5 1.0 
1.0 -0.75 0.5 1.0 


1.0 0.75 -0.5 1.0 
3.0 -0.75 -0.5 1.0 


Probfem RP4 (Configuration '. 


2 in [55]), tf = 0.25 




x < 


x > 




p u v p 


p u v p 


y>0 

y<o 


1.0 0.7276 0.0 1.0 
0.8 0.0 0.0 1.0 


0.5313 0.0 0.0 0.4 
1.0 0.0 0.7276 1.0 
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with z = x + iy, b = r^e luJt and i 2 = — 1. The circulation of each vortex is 
denoted by T, the angular rotation frequency of the vortex pair is given by 
uj = r/(47rro) and the rotation Mach number is M = r/(47rr Co), with the 
usual definition of the sound speed as Cq = JPo/Po- The ambient reference 



density and pressure are denoted by p and p 0l respectively. From (60) one 
obtains the Cartesian velocity components u x and Uy as 

dw r z 

-7- = — -2 — =v x - iv y . (61) 
dz m z 2 - b 2 y 

The hydrodynamic pressure associated with the vortex pair is given by the 
unsteady Bernoulli equation as 



V = Po - Po Re U- + - (v 2 x + v 2 v ) . (62) 




The initial density is defined by p = p 1 / 1 . Inside the vortex core radius r c , 
i.e. for r < r c , a Gaussian-type vorticity distribution is imposed according to 
[641175] . in order to avoid the singularity of the potential flow in the vortex 
center. With a matched asymptotic expansion technique [571164] . the far field 
sound pressure produced by the ideal incompressible pair of potential vortices 
can be obtained in polar coordinates r 2 = x 2 + y 2 and tan(0) = x/y as 

P' = J°Ta 2 ( J 2(^) sin (2(ut - 6)) - Y 2 (kr) cos (2(ut - 9))) , (63) 

with k = 2uj/c . 32(kr) and Y 2 (kr) are the second order Bessel functions of 
the first and second kind, respectively, and p r denotes the fluctuation of the 
sound pressure about the unperturbed ambient mean pressure p . 

For the numerical simulations, the two-dimensional computational domain of 
this problem is chosen as Q = [—500; 500] x [—500; 500] and the entire problem 
is solved with the same fourth order ADER-WENO scheme using a level zero 
grid of 250 x 250 elements, together with t = 4 and ^ max = 3. On a uniform 
fine grid this would correspond to an effective resolution of 16000 x 16000 
mesh points, hence the use of an AMR technique with local time stepping or 
at least a suitable domain decomposition with local time stepping such as the 
one presented in [98] is mandatory. For comparison, also a second order AMR 
simulation with t = 4 and ^ max = 3 is run with 500 x 500 elements on the level 
zero grid. 

The parameters used for this simulation are r c = 0.2, 7 = 1.4, p = p = 1, 



c o = y 7£>o/Po = V^> an d r = 0.08 • 47r^/7, hence the rotation Mach number 
is M = 0.08. With the above parameters the wave length of the sound waves 



A s can be computed from (63) as X s = ttcq/uj « 39. With the chosen grid 
resolution (Ax = Ay = 4 in the far field), the fourth order scheme resolves 
the acoustic waves with about 10 points per wavelength (PPW), while the 
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second order scheme employs about 20 PPW. Simulations are performed until 
t = 500, before the acoustic waves reach the corners of the outer border. The 
acoustic pressure field generated by the co-rotating vortex pair is shown in Fig. 
A comparison of our numerical simulations with the reference solution 



is depicted in Fig. 16. The reference solution (63) is a time periodic 



solution. However, it is obvious for the present problem that when starting 
from an initially undisturbed pressure field, no sound signal can arrive at a 
given spatial point x = before the time |x|/c , hence the analytical 

reference solution is depicted in Fig. 16 only for times larger than |x|/cq. Note 
further that in the present simulations the entire problem has been solved 
using the compressible Euler equations from the near field up to the very 
far field and that the singularity in the center of the potential vortices has 
been avoided by a Gaussian- type vorticity distribution inside the vortex core, 
as suggested in [6^75] . In contrast, the reference solution has been obtained 
with a matched asymptotic expansions technique for the radiated sound field 
of a pair of ideal incompressible potential vortices. Due to the modified flow 
field inside the core radius with respect to the ideal potential vortex, we expect 
our sound wave amplitudes to be always lower than the ones of the reference 
solution, which is actually confirmed by the results shown in Fig. [16] With this 
said we can observe an overall good agreement with the analytical solution 
concerning phase and amplitude of the sound pressure signal for the fourth 
order ADER-WENO scheme. For comparison, a numerical solution obtained 
with a classical second order AMR scheme on a grid refined twice as much 
is also shown in Fig. [IBj The mesh refinement for the second order scheme 
with respect to the fourth order method leads exactly to the same number 
of degrees of freedom used to represent the reconstructed solution w^. The 
second order scheme has to update four times more cells and due to the CFL 
condition also needs twice as many time steps compared to the fourth order 
scheme, which increases the number of zone updates by a factor of eight. 
Since the second order AMR scheme is 9.4 times cheaper per element update 
compared to the fourth order AMR method, see Table [TJ the total CPU times 
of both simulations are comparable. However, due to the significantly higher 
numerical diffusion of the second order method even on the refined mesh, 
an unphysical vortex merging appears, which causes the acoustic signal to 
cease completely after a certain time, since the merged vortices collapse into 
a single stationary vortex, which does not emit any sound waves. This is 
clearly seen in the acoustic signals of Fig. [TBJ which also show that the second 
order scheme obtains much lower sound pressure amplitudes in the second 
point x 2 = (200, 0) even before the unphysical vortex merging. The second 
observation point is about five propagated wavelengths away from the center 
of the vortex pair. We furthermore show the vortex configuration at the final 
time t = 500 for both the fourth and the second order scheme in Fig. 17, as 



well as the time t — 281 when the spurious numerical vortex merging takes 
place for the second order scheme. For a detailed study of physical vortex 
mergers, see [6211101 j . 
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Fig. 16. Temporal evolution of the sound pressure in the points xi = (100, 0) (left) 
and X2 = (200, 0) (right). Comparison of the second and fourth order ADER-WENO 
AMR results with the matched asymptotic expansion (MAE) solution for the far 
field sound pressure generated by an ideal incompressible co-rotating potential vor- 
tex pair. The acoustic signal of the second order scheme ceases due to a spurious 
unphysical vortex merging caused by excessive numerical diffusion. 
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Fig. 17. Density contours of the co-rotating vortex pair at time t — 500 for the 
fourth order scheme (left) and for the second order scheme (right). The spurious 
vortex merging obtained with the second order scheme is also depicted for t — 281 
(center) . 

4-2 Classical MHD equations 



In this section we consider a more complicated hyperbolic system than the 
Euler equations used in the previous section. We solve the classical, i.e. non- 
relativistic, equations of ideal magnetohydrodynamics (MHD) in three space 
dimensions. The MHD system introduces an additional difficulty for numerical 
schemes since the divergence of the magnetic field must remain zero for all 
times, i.e. 



3B X dB qi 3B Z 



0, 



(64) 



dx dy dz 

which for the continuous problem is always satisfied under the condition that 
the initial data of the magnetic field are divergence-free. From the discrete 
point of view this is not necessarily guaranteed and hence extra care is re- 
quired in the discretization. In this article we use the hyperbolic version of 
the generalized Lagrangian multiplier (GLM) divergence cleaning approach 
proposed in [27] . It consists in adding an auxiliary variable ^ and one linear 
scalar PDE to the MHD system to transport divergence errors out of the com- 
putational domain with the artificial speed c^. The augmented MHD system 
with hyperbolic GLM divergence cleaning has the state vector u given by 



u 



= (p pv T 



E B 



T 



(65) 



and the flux tensor F = (f , g, h) is defined as: 



/ 



-r 
pv 



\ 



vB - Bv + 



V 



4B T 



(66) 



/ 
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with the velocity vector v = (v x ,v y ,v z ) T ', the magnetic field vector B = 
(B x ,B y ,B z ) T and the 3x3 identity matrix I. The equation of state is the 
ideal gas law, hence 

p = ( 7 _l)(£_ lp0B_|!). (67) 



Orszag-Tang vortex system. The first test case considered for the ideal 
MHD equations is the classical vortex system of Orszag and Tang [68] which 
was studied extensively in [69] and [26J. The computational domain is ft = 
[0;27r] 2 . We use the parameters of the computation of Jiang and Wu [50] . 
scaling the magnetic field by \f^i\ due to the different normalization of the 
governing equations. The initial condition of the problem is given by 

(p, u, v,p, B x , B y ) = (V, - sin(y), sin(x), 7, -a/^tt sin(y), Vi^r sin(2a;)) , 

(68) 

with w = B z = and 7 = |. The divergence cleaning speed is set to Ch = 2.0. 
The problem is solved up to t = 5.0 using a third order ADER-WENO scheme 
with componentwise WENO reconstruction. The initial mesh on level zero is 
composed of 50 x 50 elements. We furthermore use t = 4 and ^ max = 2. This 
corresponds to an equivalent resolution on a uniform fine mesh of 800 x 800. To 
assess the accuracy and efficiency of our proposed AMR scheme, we also run a 
simulation on the uniform fine mesh for comparison. The results for pressure 



are shown in Fig. [18] for t = 0.5, t = 2.0, t = 3.0 and t = 5.0, both, for the 
AMR grid as well as for the uniform grid corresponding to the finest AMR grid 
level. Our results are in agreement with the fifth order WENO finite difference 
solution computed by Jiang and Wu [50] and with the unstructured third order 
WENO solution depicted in [28] for the same problem. Furthermore, the AMR 
computations are in excellent agreement with the uniform fine grid reference 
solution. An efficiency comparison concerning memory requirements and CPU 
time is listed in Table [5j The AMR method needs 322 time steps on the coarsest 
mesh, while the simulation on the uniform fine mesh needs 4148 time steps 
to reach the same final time. Even for this test problem, where most of the 
cells are refined and therefore only little gain is expected through the use of 
AMR techniques, we obtain still a speedup of a factor of 1.8 compared to the 
uniform fine mesh simulation. 

Table 5 

Memory and CPU time comparison of the third order ADER-WENO AMR method 
and ADER-WENO on a uniform fine grid for the Orszag-Tang problem. Memory 
consumption is measured in maximum number of elements and CPU time is nor- 
malized with respect to the simulation on the fine uniform mesh. 





AMR 


Uniform 


ratio 


Cells 


454525 


640000 


1.41 


CPU 


0.547 


1.0 


1.83 
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Fig. 18. Orszag-Tang vortex system at times t = 0.5, t = 2.0, t = 3.0 and t = 5.0 
from top to bottom. AMR grid (left), third order ADER-WENO solution obtained 
on the AMR grid (center) and on a fine uniform grid corresponding to the finest 
AMR grid level (right). 

MHD rotor problem. The second test case is the well-known MHD rotor 
problem proposed by Balsara and Spicer in [7] . It consists of a rapidly rotating 
fluid of high density embedded in a fluid at rest with low density. Both fluids 
are subject to an initially constant magnetic field. The rotor causes torsional 
Alfven waves to be launched into the fluid at rest. As a result the angular 
momentum of the rotor is diminished. The problem is set up on a computa- 
tional domain Q = [-0.6; 0.6] x [-0.6; 0.6], using a third order ADER-WENO 
scheme. The AMR mesh on level zero contains 60 x 60 elements. With t = 4 
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and ^ max = 2 this simulation corresponds to a uniform fine mesh with 960 x 960 
points resolution. As before, a uniform fine grid simulation is also performed 
to assess accuracy and efficiency of the proposed ADER-WENO scheme on 
AMR grids. The initial density of the rotor is p = 10 for < r < 0.1 and 
p = 1 for the ambient fluid. The rotor has a constant angular velocity u that 
is determined in such a way to obtain a toroidal velocity of v = uj • r = 1 at 
r = 0.1. The pressure is p = 1 in the whole domain and the magnetic field 
vector is set to B = (2.5, 0, 0) T in the entire domain. As proposed by Balsara 
and Spicer we apply a linear taper to the velocity and density field, however 
only in a very small range 0.1 < r < 0.105 so that density and velocity match 
those of the ambient fluid at rest at a radius of r = 0.105. The speed for the 
hyperbolic divergence cleaning is set to = 2 and 7 = 1.4 is used. Transmis- 
sive boundary conditions are applied at the outer boundaries. The final AMR 
mesh is depicted in Fig. [20) The computational results on the AMR mesh are 



compared with those on the uniform fine mesh at time t = 0.25 in Fig. [19 
for density, pressure, Mach number and magnetic pressure. One observes that 
both solutions agree very well with each other. Also compared to the results 
presented by Balsara and Spicer we note a very good agreement. We emphasize 
that thanks to the divergence cleaning, no spurious oscillations can be seen 
in the density field and in the magnetic pressure, as reported by Balsara and 
Spicer for Godunov schemes without divergence cleaning. The AMR method 
needs only 99 time steps on the coarsest mesh, while the simulation on the 
uniform fine mesh solution needs 1147 time steps to reach the same final time. 

In this test problem, the efficiency gain of AMR is particularly evident. The 
computation on a fine uniform mesh corresponding to the finest AMR level 
needs more than five times more elements and more than seven times more 
CPU time, see the detailed results reported in Table [6) 



5 Conclusions 



In this article we have presented the first better than second order one-step 
ADER-WENO finite volume scheme on space-time adaptive AMR grids. The 
use of a high order one-step time stepping method, based here on a local 
space-time discontinuous Galerkin predictor, allows a straightforward imple- 

Table 6 

Memory and CPU time comparison of the third order ADER-WENO AMR method 
and ADER-WENO on a uniform fine grid for the MHD rotor problem. Memory con- 
sumption is measured in maximum number of elements and CPU time is normalized 
with respect to the total wallclock time on the uniform mesh. 





AMR 


Uniform 


ratio 


Elements 


179680 


921600 


5.13 


CPU time 


0.140 


1.0 


7.14 
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Fig. 19. MHD rotor problem at time t = 0.25. Third order ADER-WENO solution 
obtained on the AMR grid (left) and on a fine uniform grid corresponding to the 
finest AMR grid level (right). 
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Fig. 20. AMR grid for the MHD rotor problem at time t = 0.25. 



mentation of time accurate local time stepping, where each AMR grid level 
runs on its own local time step. Furthermore, compared to the method of 
lines based on Runge-Kutta time stepping, the use of a high order one-step 
scheme in time reduces the number of nonlinear WENO reconstructions and 
the number of necessary MPI communications. All these key features of our 
present scheme help to keep the overall overhead associated with the admin- 
istration of the space-time adaptive mesh at a reasonable level, at most 25%, 
as quantified in Table [TJ 

We have carried out numerical convergence studies, confirming that the claimed 
higher order in space and time is actually reached in practice. Furthermore, 
the scheme has been applied to a series of test problems in two and three space 
dimensions, solving the compressible Euler equations as well as the classical 
MHD equations. In our examples, it was clearly shown that also for better 
than second order schemes, the use of AMR is beneficial, compared to the use 
of a uniform fine grid. We have also shown via numerical evidence, that even 
in the AMR context the use of higher order schemes is beneficial, in particular 
when small scale turbulent structures, vortices and sound waves have to be 
resolved. All these physical phenomena require little numerical dissipation for 
their efficient simulation. 

Future research will concern the extension of the present method to the general 
family of PnPm schemes introduced in [28] . as well as the simulation of realistic 
problems in computational astrophysics. In further work we plan to extend 
the present scheme also to turbulent viscous flows, to chemically reacting 
multiphase flows as well as to nonconservative hyperbolic systems. 
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